Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 2 November 20 1 2 (MN KT^C style file v2.2) 



The Effects of Radial Migration on the Vertical Structure of 
Galactic Discs 

Rok Roskar 1 *, Victor R Debattista 2 , Sarah R. Loebman 3 

1 Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland 

2 Jeremiah Horrocks Institute, University of Central Lancashire, Preston, PR1 2HE, UK 

3 Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA 



12 November 2012 



ABSTRACT 



We present evidence that isolated growing discs, subject to internal spiral perturbations, 
thicken due to both heating and radial migration. We show this by demonstrating that the 
thickness and vertical velocity dispersions of coeval stars depend on their age as well as the 
change in their radii. While the disc thickens due to internal processes, we find that this in- 
duces only a minor amount of flaring. We further demonstrate the consequences of such thick- 
ening on the structural properties of stellar populations and find that they qualitatively agree 
with recent studies of the Milky Way disc. 



1 INTRODUCTION 



Stars born at one radius within a disc galaxy such as the Milky 
Way (MW) need not remain at that radius indefinitely, even if they 
retain nearly circular orbits. Sellwood & Binney (2002 ) have shown 
that the interaction of stars with transient spirals at their corotation 
resonance can lead to substantial changes in angular momentum of 
the stars without significant heating. This radial migration has many 
important consequences for properties of stellar populations, such 
as giving rise to a flattened age-metallicity relation with substantial 
scatter ([Sellwood & Binney|2002[ [Roskar et al.|2008a||Schonrich| 
|& Binney|2009a| ). 

While it is typically assumed that the vertical and in-plane mo- 
tions of stars are largely decoupled, it is becoming increasingly ap- 
parent that some coupling exists and affects the detailed chemody- 
namical structure of discs. For example, as stars migrate outwards 
they experience a smaller vertical restoring force and are therefore 
expected to reach greater height above the mid-plane. Schonrich & 
Binney (2009a b) explored this idea by using analytical chemical 
evolution models that allowed for stars to radially migrate and as- 
sumed that in the process of migration their vertical energies were 
conserved. They showed that under these assumptions it was pos- 
sible to create a thick disc entirely from migrated stars that shared 
many of the chemodynamical properties of the thick disc in the 
MW. 

Disc thickening due to migrating stars has also been studied 
in simulations by |Loebman et al.| ( [201 1| ), who showed that out- 
wards migrating stars end up in the more vert ically extended com 



ponent of a double seen profile density fit. As Schonrich & Binney 



(2009b|, they found that cuts by kinematics or abundances gives 
rise to an artificial thin/thick disc dichotomy. Thus, radial migration 
appears as yet another possible mechanism for creating a thickened 
disc component, but because it relies entirely on internal processes 
it is conceptually orthogonal to most other thick disc formation sce- 
narios. 



Other mechanisms for creating a thickened disc component 
involve either heating the disc through bombardment (e.g. Quinn 
|et al.|1993] >, forming a thickened disc from the remnants of a gas- 
rich last major merger ( [Brook et al.pOO^ or accreting an extra- 



planar stellar component ( [Abadi et al.| |2003). All three of these 



mechanisms are a result of the cosmological environment and con- 
sequently the thick disc is often considered to be an indicator of the 
MW's cosmic history. On the other hand, a thickened component 
built up by migration has little to do with cosmological evolution 
but depends largely on internal disc dynamics. If migration can lead 
to thickening of the stellar component then it can muddle the sig- 
natures of cosmologically-relevant events that undoubtedly marked 
the early stages of the MW's disc formation. The formation of the 
thick disc through heating by clump instabilities is in principle also 
an in-situ formation mechanism (Bournaud et al. 2009 ), but it still 
relies on the cosmological environment for the rapid gas accretion. 

The suggestions that the MW thick disc may partially form 
through migration of stars forming in the thin disc led to a num- 
ber of observational studies exploring the formation of the thick 
disc. Sale s et ah] ([2009 ) proposed using the eccentricity distribu- 
tion of thick disc stars to distinguish between thick disc forma- 
tion mechanisms. Consequently, several groups sought to apply the 
Sales et al. ( 20 09} eccentricity test to observa tional data. jDie rickx 
|et al.| ( |2010| ), |Casetti-Dinescu et al.| ( [20TT] ), and |Wilson et al.| ( |201 1| ), 
all reached a general conclusion that the eccentricity distribution is 
mostly inconsistent with the accretion scenario, but broadly con- 
sistent with the other three (the gas-rich merger scenario being the 
most favorable). The eccentricity distribution is not alone enough 
to break the degeneracies among the models, but such constraints 
become stronger once metallicities and abundances are also con- 
sidered: |Liu & van de Ven] ( |2012|) showe d that low eccentricity 
stars from a SEGUE survey ( [Yanny et al.||2009| > sample follow a 
continuous distribution in the [Fe/H]-[a7Fe] plane while the high 
eccentricity stars seem disconnected in this plane. In contrast Bovy 
|et al.| ( [201 2a| ) found little evidence of a distinct thick disc in ei- 
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ther a mass-, or number-, weighted distribution of stars from the 
SDSS/SEGUE dataset. 

Separated into mono-abundance populations, a-rich, metal- 
poor, large scale-height populations in the MW have a shorter 
scale-length, while metal-rich stars have low scale-heights and long 
scale-lengths ( |Bovy et al.|2012c|> (a similar result was found for the 
high-a population by |Cheng et al.||2012| ). The stars in individual 
abundance bins are also vertically isothermal (Bovy et al. 2012b ), 
strengthening the argument that they are single populations. Both 
of these properties seem to favour strong internal evolution for the 
MW disc, because if a disc builds up entirely from internal pro- 
cesses, a distinct second component that dominates away from the 
plane does not readily form. It is, however, also possible to recover 
many of these trends in the cosmological context if the disc ini- 
tially forms hot and thick, subsequently forming increasingly thin- 
ner populations as its self gravity gradually increases (Stinson et al. 
in prep). 

However, even if stars form in a thicker component early in 
the history of the disc, they do not become immune to disc per- 
turbations. iV-body simulations have shown that stars in the thick 
disc were also able to migrate though to a lesser extent (Solway 
|et al.|2012] >. The same study showed that migrating stars on aver- 
age (but not individually) conserve their vertical action, rather than 
their vertical energy. 

Recently, Schonric h & Binney| ( |2012| > studied the co- 
dependence of vertical and horizontal motions analytically. They 
found that the vertical motion significantly affects the detailed in- 
plane velocity distributions. Following Binney ( 2010 ), the coupling 
between the two components of motion was achieved through the 
assumption of vertical action invariance. This assumption has been 
shown to be a reasonable approximation by more rigorous torus 
modelling (Binney & McMillan 2011). While the above analytic 
studies have focused on the co-dependence of oscillations about 
the guiding centre and the vertical motion, in this Paper we explore 
the implications of radial migration (i.e. the change in the guiding 
centres) for the vertical distribution. Under action conservation in a 
radially-decreasing disc potential, any change in radius will lead to 
a modification of the vertical motion. However, while stars oscil- 
late in radius by at most up to a kpc they might migrate radially by 
many kpc and we can therefore expect that migration would have a 
larger effect on their vertical motion. 

However, Min chev et al.| ( |2012a| ) argue that vertical action 
conservation prevents the stars from thickening as they migrate 
outwards. They support this claim with results from isolated N- 
body+SPH simulations from the GalMer database ( Chilingarian 
|et al.|20T0|> as well as s ticky-particle semi-cosmological simula- 
tions Martig et al. 2012|). Their conclu sion disagrees sharply with 
the findings of Schonrich & Binney (2009b), who assumed en- 
ergy conservation rather than action conservation when modeling 
the vertical distribution in their analytic models. However, it also 
disagrees with Loebman et al.| ( |2011| ) who used high-resolution 
iV-body/SPH simulations and found that the migrated population 
formed a thicker component that was shown to be broadly consis- 
tent with MW thick-disc trends in both chemistry and kinematics. 

In |Roskar et al.| ( [2bl2| ), we focused on the details of the ra- 
dial migration mechanism in a suite of iV-body/SPH simulations. 
We found evidence that the largest migrations were taking place at 
the corotation resonance (CR) of dominant spirals, confirming the 
mechanism proposed by Sellwood & Binney ( 2002 ). In this paper, 
we elucidate the connection between radial migration and vertical 
disc structure using the fiducial simulation of Roskar et al. ( 2012| ). 
To gain insight into the physical nature of the processes involved, 



we compare the simulation data with expectations based on the sim- 
plified analytic models of Schonrich & Binney (2012 ). The paper 
is organized as follows: our methods are presented in Section [2] 
the basic results are presented in Section]!] we compare the results 
to a simple analytic model for disc thickening in Section]!] we dis- 
cuss the consequences of such thickening on the trends in structural 
properties of stellar populations in Section[5] 



2 SIMULATIONS 

We analyze the same iV-body/SPH simulation that we used in 
Loebman et al. ( 2011}. This simulation is a re-run of the fiducial 
simulation analyzed in |Roskar et a l. (2008b a 2012 ) using metal 
diffusion, but is otherwise identical in every respect. The initial 
conditions consist of two NFW ( Navarro et al. 1997 ) halos, one of 
dark matter with a mass of 10 12 M and the other of gas in hydro- 
static equilibrium with a mass of 10 11 M©, both sampled with 10 6 
particles. The gas also rotates with A = 0.039 and j oc R, where 
j is the specific angular momentum and R is the cylindrical radius. 
We use the code GASOLINE ( Wad sley et al.|2004| ) to perform the 
computation. Once the simulation begins, the gas cools and col- 
lapses into a disc that forms stars according to a standard prescrip- 
tion ( [Stinson et al. 2006 ) with a density threshold of 0.1 amu/cc and 
a temperature cutoff of 1.5 x 10 4 K. Star particles provide "feed- 
back" to the gas by injecting it with energy and polluting it with 
supernova ejecta. The simulations are evolved completely in iso- 
lation and no cosmological context is included. By the end of the 
simulation, the disc is populated by more than 2 x 10 6 star par- 
ticles. The softening we use for the baryonic component is 50 pc. 
The reader is referred to Roskar et al. (2008b ) and Roskar et al. 
( 2012 ) for further details on the simulations. The simulations were 
analyzed with the aid of the community-developed PYTHON-based 
analysis package Pynbod^ 

Our simulation code includes prescriptions for the generation 
of metals in supernova type la and II explosions as well as in AGB 
stars. We can therefore follow the relative abundances of a ele- 
ments relative to Fe (Oxygen is used as the a element proxy). Due 
to uncertainties in metal yields, the absolute values of [O/Fe] do 
not match with observations, but relative trends are nevertheless in- 
formative. For details on the metal enrichment implementation see 
Stins on et al.| ([2006). As in Loeb man et al.| ( [20TT| ), the simulation 
we use here also implements a prescription for the diffusion of met- 
als through the ISM. 



3 RESULTS 

The simulation yields a disc that is kinematically similar to that of 
the Milky Way. In Fig. [T] we show the age-velocity dispersion re- 
lation for stars that are found in the "solar neighborhood" at the 
end of the simulation. We define solar neighborhood as a disc re- 
gion between 7.5-8.5 kpc and within 200 pc from the plane. The 
disc scale length in our model is ~ 2.5 kpc, so this is a reason- 
able approximation to the region of the MW disc where analogous 
relations have been observed. The age-velocity relation from the 
Geneva-Copenhagen survey (Holmberg et al. 2009 ) shows similar 
power-law dependencies and values, though in our case the oldest 
stars are hotter by ^ 10%. Part of this increased dispersion may be 
due to the fact that our simulation yields a slightly more massive 
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Figure 1. Age-velocity dispersion plot for stars with 7.5 < R [kpc] < 8.5 
and \z\ < 200 [pc] after 10 Gyr of evolution. 

disc (the circular velocity v c ~ 250 km/s). Fitting power-laws to 
the relations in Fig. [T] we obtain power-law indices of 0.24, 0.26, 

0. 25, and 0.15 for a tot, ctr, cr<f>, and a z respectively. These indices 
are lower than those observed in the MW, indicating that an addi- 
tional source of heating may be needed. On the other hand, the fact 
that the heating indices are not too high also means that the am- 
plitudes of non-axisymmetric structure forming in the disc are not 
unreasonably high. A related issue is that the youngest stars are not 
quite as kinematically cool as the observations due to the dispersion 
floor in the gas component, which also contributes to higher overall 
velocity dispersions. We discuss this further in Sect. [4] 

3.1 Thickening and vertical velocity dispersion due to 
migration and heating 

Because radially migrating stars feel a smaller restoring force to- 
wards the midplane as they move outward, the amplitude of their 
vertical oscillations increases. However, at the same time, heating 
processes due to the recurring spiral structure also heats the stars, 
regardless of whether they have migrated or not, providing another 
source of increased vertical motion. Comparison between Figs. 3 
and 5 in Loeb man et al.|p0TT) shows this implicitly. Their Fig. 3 
shows that when significant migration is present, the stars dominat- 
ing the distribution away from the plane have also migrated from 
the interior of the disc. By contrast, in the model with less migration 
shown in their Fig. 5, the populations away from the plane are there 
mostly because of heating (i.e. they are old) rather than migration 
(they formed in- situ). 

In Fig. [2] we show the relative dependence of stellar popu- 
lation thickness on radial migration and heating. As a measure 
of thickness, we use the model-independent z r ms = Each 
panel shows the distribution of Az rms = z r ms,now — z rms ,form, 

1. e. the change in thickness since birth, as a function of age and 
AR = Rnow — Rform, for a given range of formation radii. The 
distributions shown in each of the panels are made by selecting par- 
ticles in a particular range of formation radii (indicated at the top of 
each panel). We bin the particles in this space on a grid - each bin 
then corresponds to roughly a coeval population (born at approxi- 
mately the same radius, of the same age, and migrated by the same 
amount). The contours show the mass density of particles. 

If all the changes to the thickness of each population were due 
to heating, we would expect age to be the determining factor in set- 



ting the Zrms- As a result, the 2D gradient in Az r ms would be pre- 
dominantly in the vertical direction. Conversely, if all the changes 
in thickness were due to radial migration, the gradient would be 
horizontal. 

The gradients in all of the panels of Fig. [2] are diagonal, with 
Azrms increasing in the direction of older ages and larger AR, 
implying that radial migration as well as heating contribute sub- 
stantially to the vertical thickness of the stellar populations. This is 
true for all ages and all AR. For positive AR, Az r ms is a mono- 
tonic function of AR. Stars migrating inward decrease their z rms 
so long as they stay away from the central bulge region. Those stars 
that migrate to the very centre of the disc increase their ^ ms , pre- 
sumably due to rapid heating that occurs there due to the presence 
of multiple inner Lindblad resonances and a weak oval that devel- 
ops in the centre from time to time. The high velocity dispersions 
for these stars indicate that this may indeed be the case (see Fig. [3] 
discussed below). The oldest stars in the interior of the disc (upper- 
left corners of the top two panels) are on eccentric orbits and there- 
fore not as affected by disc perturbations but their z r ms decreases 
due to adiabatic contraction. Our main finding from Fig. [2] is that 

1) at any given age, thickness is a monotonic function of AR, and 

2) at any AR thickness is a monotonic function of age. 

We find a very similar kind of co-dependence for a z , the ver- 
tical velocity dispersion, shown in Fig. [3] The velocity dispersions 
of the stars migrating outward are lower than the dispersions of 
those migrating inwards. This is to be expected because the stars are 
moving to a region of a decreased midplane potential and because 
stars in the inner disc heat very efficiently (due to inner Lindblad 
resonances or non-axisymmetric structure). Nevertheless, the dis- 
persion clearly depends on both parameters, age as well as distance 
migrated since birth, consistent with Fig. [2] 

In Figs. [2] and [3] we have shown that the thickness and verti- 
cal velocity dispersion of coeval stellar populations depend on the 
magnitude of migration as well as their age. In other words, at a 
fixed present-day radius, the velocity dispersion and thickness of a 
population depends on its birth location and its time of formation. 
However, we have also shown that a z of stars actually decreases as 
they migrate outward, so can these stars entering the solar neigh- 
borhood from the inner disc then still masquarade as the thick disc? 
From Fig. [3] we can see that the old (> 8 Gyr) stars that have mi- 
grated outward by several kpc have vertical velocity dispersions of 
40-50 km/s. These stars comprise the a-old population that in the 
MW has very similar velocity dispersions ( |Liu & van de Ven|2 012 ; 
|Bovy et al.|2012b) ; note that for the purposes of defining the thick 
disc, the "characteristic" a z is taken to be ~ 35 km/s ( Bensby et al. 
2003), a criterion met by this migrated population. Even if the pop- 
ulations "cool" as they migrate outwards, their present-day prop- 
erties can remain consistent with a thickened, a-rich population. 
Note that in our models the only stars with the appropriate a z are 
> 5 Gyr old and certainly the only ones that are thickened enough 
to represent the thick disc are even older and have migrated sub- 
stantially. The in- situ young stars have a z < 20 km/s, making up 
the "thin" disc. 

Minc hev et al.| ( |2012a| ) similarly found that velocity disper- 
sions decrease as stars migrate outwards. However, they interpreted 
the decrease in a z as an indication that radial migration cannot con- 
tribute to a thickened population. We find instead that it is indeed 
possible to form a thickened and kinematically hotter component 
with the aid of radial migration. However, because the stars migrat- 
ing from the inner to outer disc begin with a relatively high veloc- 
ity dispersion, they still end up kinematically hotter relative to the 
in-situ population when they arrive to the outer disc. Furthermore, 
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AR [kpc] AR [kpc] 

Figure 2. Mean change in thickness as a function of age and AR plotted in four different bins of formation radii. Contours show the particle mass density 
and are logarithmicaly spaced from 10 to 10 4 particles per bin. Colours correspond to the change in thickness, defined as the root-mean-square of the vertical 
displacement from the midplane, z rrns = V¥*. Increase in thickness is a function of both, age and AR, i.e. stars undergo heating but their vertical distribution 
is also affected by the effective change in potential as they migrate inwards or outwards in the disc. Note that the x-axis range is different in each panel. 



as Fig. [2] clearly shows, they thicken vertically as well. We spec- 
ulate the discrepancy between our findings and those of Mi nchev| 
et al. (2012a) is partly due to the rapid heating apparent in their 
simulations. Shortly after their runs are started velocity dispersions 
increase dramatically, indicating that those models are initially vi- 
olently unstable. Beyond this initial instability, the phenomenon 
driving the disc trends in their models (the same simulations were 
also discussed in Minchev et al. 2011 and Minchev et al. 2012b ) 
is rapid angular momentum redistribution during bar growth (e.g. 
Hohl 1971 ) as opposed to redistribution by spirals which dominates 
in our simulations ( [Roskar et al.|2012j >. 



3.2 Disc flaring 

An important consideration for any model that deals with the thick- 
ening of the stellar disc is the change in scale height with radius, 
i.e. disc flaring. While substructure bombardment is particularly ef- 
ficient at thickening discs, it may also cause rather drastic flaring to 
occur (Kazantzidis et al. 2009 ). Such flaring is not typically seen 
in extragalactic observations of edge-on discs ( de Grijs & van der| 
|Kruit|1996| ). 

While, radial migration has also previously been assumed to 
cause dramatic flaring (Minch ev et al.|20 12a), we show in Fig. [4] 
that this is not necessarily the case. In the top panel, we show the 
vertical profiles of the disc at several different radii, all normalized 
to their values at 1 kpc from the plane to more easily compare the 
profiles furthest from the plane. The vertical distribution clearly 
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Figure 3. Same as Fig.[2]except here the colours correspond to a z , the vertical velocity dispersion at the end of the simulation. The old stars that the "solar 
neighborhood" (i.e. ~ 8 kpc) from the interior of the disc have a z ~ 40 — 50 km/s in broad agreement with thick disc values in the MW. 



becomes more extended in the outer parts of the disc, but the flaring 
is rather minor. 

In the bottom panel, we quantify the flare by showing the val- 
ues of maximum-likelihood estimates (MLE) for the scale height of 
the thicker component of a double sech 2 distribution as a function 
of radius. We use MLE instead of a least- squares fit to the binned 
vertical profiles because we found those fits to be unreliable, while 
MLE gives consistently good results. To find the MLE parameters, 
we maximize the log-likelihood function 

N 

l(h 1 ,h 2 ,f\z) = J2^[p(Mhi,h 2 J)] - Hpo), 



with 



p(z) = (1 - /) sech 2 0/2/n) + / sech 2 0/2/i 2 ) 



where p is the volume density, N is the total number of particles 



in the bin, z% is the vertical position of the i-th particle, hi and 
/12 are the thin and thick disc scale heights, / is the thick disc 
fraction, and po is the midplane density. To estimate the parame- 
ters, we select particles found in 1 kpc radial bins centered on each 
point shown in Fig.[4]and up to 4 kpc from the midplane. Following 
|Bovy et aL| ( |2012c|>, the likelihood function I is maximized using an 
MCMC sampler ( Foreman-Mackey et al. 2012), yielding parameter 
estimates and uncertainties, represented by the error bars in Fig. [4] 

The increase in scale height is only 50% across approximately 
10 kpc, corresponding to 4 disc scale lengths. The thick disc frac- 
tion is rsj 12 — 20% in the main disc region (4-10 kpc) and rises 
to ~ 40% in the outer disc. Interior to 10 kpc, the scale height 
is essentially constant. Note that most studies fitting the vertical 
density of the Milky Way disc (e.g. |Gilmore & Reid||1983| |Juric| 
|et al.|2008||Bovy et al.|20T2c] use an exponential fitting function 
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Figure 4. Top: Vertical density profiles at several different radii, normalized 
to their values at z = 1. Bottom: Scale heights of the seen 2 component as 
a function of radius. From the inner disc to the outermost regions we find a 
modest increase in scale height. 



and that the scale heights between the two functional forms are not 
strictly comparable. We find that fitting our simulated particle dis- 
tribution with an exponential model yields scale heights larger by 
up to ~ 30%, making them only slightly smaller than similar thick 
disc scale height measurements in the MW. However, we use the 
sech 2 because it gives a better fit to our model. 



4 DISC THICKENING DUE TO CONSERVATION OF 
VERTICAL ACTION 

Here we consider the simplified ID problem with a star oscillat- 
ing above and below the plane of a thin disc. This is similar to the 
model described by Minc hev et al.| ( |2012a| ), who use it to show that 
the vertical velocity dispersion is expected to decrease as stars mi- 
grate outward if the vertical action is conserved. We show that this 
same model leads to the conclusion that the stellar vertical distri- 
bution must thicken. In this simple case, the vertical frequency is 



v = 2 v / ttG'£0R), (1) 

where S(i?) is the surface density of the disc at radius R and G is 
the gravitational constant. The vertical action can be approximated 



by L z — E z /v, where E z is the vertical energy. For a population 
of stars we assume th at E z ~ a 2 . If we assume that the vertical 
action L z is conserved { [Solway et al.|2012| then a 2 oc v. Assuming 
oc e~ R/Rd , where Rd is the scale length, Eqn.[l]yields 

a 2 oce- R/4R \ (2) 



as shown previusly by Minchev et al. (2012a). 

The vertical density distribution of an isothermal population 
is given by ( Spitzer 1942 ) 

p(z) = po sech 2 (z/2h z ), 



with the vertical scale-height 



(3) 



(4) 



* 2ttG£ 

Plugging Eq.^into Eq. [4] one finds that h z oc e R/2Rd , i.e. the 
scale height increases with positive changes in radius. Combining 
the approximation for change in a z for a migrating population with 
Eq.[4]yields a simple expression relating the change in radius AR 
and the ratio of initial and final scale heights, h z j/h z ,i, and Rd'. 



hz,f 
h z .i 



exp 



(AR\ 



(5) 



\2RdJ 

This simple derivation of Eq. [5] ignores the presence of other disc 
components; such an isothermal population is just one of many that 
make up the disc. In particular, the addition of mass (by gas accre- 
tion and star formation) in the midplane would adiabatically com- 
press the migrated population, which has previously been shown to 
alleviate some of the thickening effects of infalling substructures 
(Villalobos et al. 2010 ). Furthermore, even in a disc not subject to 
outside perturbations, <j z will also increase with time due to heat- 
ing from disc structure thereby further boosting the scaleheight, an 
additional effect that we discuss below. 

Eq.|5]is an approximation to the more general expression de- 
rived by Schonrich & Binney (2012) (their Eq.[27]), which instead 
of the factor 2 on the right hand side includes a factor of (2 + a). 
The a changes depending on the magnitude of vertical oscillation. 
Note that the change induced in the z distribution due to adiabatic 
invariance occurs irrespective of the reason for the change in radius. 
The change may be due to radial oscillations about the guiding ra- 
dius ( Schonrich & Binney 2012 ) or due to the change in the guiding 
radius itself. Therefore, the increase in z r ms with increasing AR in 
Fig.[2]is due to a complex combination of radial heating, migration, 
and the evolution of the galactic potential. 

In Fig. [5] we try to separate out the dependence of stellar pop- 
ulation thickness on these processes. We focus on stars that are 5-6 
Gyr old and formed within the range of specified radii. We then 
further separate these stars by their orbital circularity parameter, 
j z /j c (E) measured at the present time, where j z is the specific 
angular momentum and j c (E) is the maximum specific angular 
momentum for a star with a given energy. The solid lines show 
the variation in fractional change in thickness z r ms/ 'zi^ms versus 
change in guiding radius since birth, AR g for different ranges in 
orbital circularity. The range in j z /j c is quite representative of this 
particular subset of stars, i.e. there are not many particles within 
the chosen constraints that are significantly kinematically hotter. 
We define R g in terms of angular momentum, i.e. R g = j z /V c , 
where V c = 250 km/s is the approximate circular velocity. The 
age and Rform ranges are chosen to give a reasonable sampling of 
the main part of the disc. 

Looking at Fig. [5] one can now identify the dependence of a 
stellar population's vertical distribution on heating and migration 
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Figure 5. The ratio of z r ms (thickness now) to Zi, rrns (thickness at formation) of 5-6 Gyr old stars that were formed within the indicated radial range, plotted 
against the change in guiding radius AR g = Aj z /V c from their birth to the present, where Aj z is the change in angular momentum and V c is the circular 
velocity. The colours of the solid lines correspond to different ranges in the orbital circularity parameter j z /j c (E). The overplotted red dashed lines are 
exponential fits, indicating that the changes in thickness satisfy the expectations from adiabatic invariance of the vertical action (see text). The other coloured 
dashed lines show the relation from Eq.[27] of Schonrich & Binney ( 2012) for the kinematically coldest population, overplotted for different values of a and 
assuming a scale length of 2.5 kpc. 



separately. For the kinematically coldest population, the thickness 
changes by a factor of 3 along the range of AR g . On the other hand, 
it increases by a further factor of ~ 2 from the coldest population 
to the hottest population, indicating that both processes actually 
contribute in similar amounts to the thickening. 

The evolution of the thickness of this coeval stellar popula- 
tion can be described by the simple disc model discussed above 
under the assumptions of adiabatic invariance. We extend the sim- 
ple formula in Eq.[5]to that derived by Schonrich & Binney ( 2012 ) 
and plot the result in Figure [5] for different values of a shown with 
coloured dashed lines using a disc scale length Rd = 2.5 kpc. The 
least- squares best fit is shown with the red dashed line. Note that 
for a sech 2 distribution, z r ms is smaller than the scale height by 
approximately 10%, so z r ms gives a reasonable approximation to 
zo that we can use in Eq. [5] 

The agreement between the simulation and the analytic pre- 
diction of the simple model, given the assumptions, is encourag- 
ing. Note that the likely value of a should be 1-2 since particles 
with the largest values of j z /jc have amplitudes of vertical oscilla- 
tions \zmax | <C Rd (Schonrich & Binney 2012), which provides an 
excellent description of the simulation data. We have experimented 
with using different criteria for the Rform and age selection and 
found no appreciable differences, though the best- fit parameters 
start to deviate from the expected values close to the centre or far 
in the outer disc where the surface density profile can no longer 
be described as exponential. The three different ranges of Rform 
shown in Figure[5]show no appreciable differences. If we choose a 
different age range, the relations stay qualitatively the same. 

Heating by spirals in our simulation results in increased ran- 
dom motion in-plane and out-of-plane. The increased thickness of 
the populations on slightly more eccentric orbits shown with the 
blue and green lines is therefore due to the combined effects of 
higher vertical random energy as well as amplified vertical oscilla- 
tions due to adiabatic invariance at the orbital apocentres. 

We have shown above that the thickness of a coeval popula- 
tion is an exponential function of its migration in the disc. It is 
interesting to pose a similar question for stars that share the same 
angular momentum today but that may have originated in differ- 
ent parts of the disc, i.e. are not to be considered a coeval popula- 



tion. As in Fig. [5] in the leftmost panel of Fig. [6] we show the ratio 

Zrms/zi,rms for stars with 5 < Rfinai [kpc] < 5.5 vs. AR g . 
We see similar exponential trends at all values of j z /jc, as before. 
Since the most circular (red line) population is the least affected by 
heating, we can once again use it as a proxy for the overall effect 
of radial migration on population thickness. Across the range of 
ARg, we find that migration can lead to an thickness increase by a 
factor of 2.5 (we cannot eliminate all heating so thickness increases 
by 50% simply by virtue of these stars spending 5 Gyr in the disc). 
Again, taking the range from blue to red as the effect of heating, we 
see that heating contributes another factor of < 1.5. Based on this 
and Fig. [5] we conclude that migration thickens stellar populations 
by about the same proportion as the internal heating processes, for 
parts of the disc outside the bulge. 

In our simulations, the gas cooling function does not include 
metal line or molecular cooling, limiting the temperature of the gas 
to > 10 4 K. As a result, the gas velocity dispersion is largely set 
by the minimum temparature and consequently the gas scale height 
increases as a function of radius. Due to this gas flaring, the popu- 
lations at each AR 9 in Fig.[6]are actually born with different thick- 
nesses, as shown in the middle panel. The variation across the range 
of 5-6 Gyr old stars that end up at the specified angular momen- 
tum is a factor ^2 — 3. Hence, while it is inevitable that individ- 
ual populations change their vertical distribution as they migrate 
(Fig. [2]), these changes are partially offset by the flaring of stars at 
birth (middle panel of Fig. |SJ. As a result, at a fixed final angular 
momentum (i.e. R g ) there is little variation of thickness with AR g . 

However, the molecular and atomic gas components in the 
MW show little variation in thickness out to ~ 10 kpc ( |Bronf-| 
|man et al.|1988l|Wouterloot et al.|1990[|Narayan & Jog|2002| ). We 
may therefore expect that in the MW the flaring would not mask 
the thickening resulting from radial migration. Evidence for simi- 
larly constant scale heights as a function of radius has also recently 
been found with resolved-star studies of edge-on nearby discs (R. 
de Jong & D. Streich, private communication). 

We checked the results of this section against a higher- 
resolution version of our simulation that uses four times as many 
particles (simulation R4 from Ros kar et al.| |2012) and found no 
significant differences, so we consider the dynamical results to 
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Figure 6. The change in population thickness as a function of AR 9 , the change in guiding radius since birth, defined as R 9 = j z /V c (same as in Figure[5j 
for particles with R g in the indicated range. The left panel shows the ratio of z r ms / Zi, rrns , predicted by Eq.[5]to follow an exponential. The centre and right 
panels show z r ms at time of formation and at the present time respectively. Stars used in this figure are 5-6 Gyr old. 



be numerically robust. However, assessing the significance of the 
gas disc flaring and the subsequent radial dependence of the scale 
height of young stars is a much more difficult problem. It depends 
sensitively on the details of the gas physics on scales that are an or- 
der of magnitude smaller than what we (or any other similar state- 
of-the-art simulations at present) are able to adequately resolve. 



5 CONSEQUENCES FOR LOCAL ABUNDANCE 
TRENDS 

The consequence of the thickening described above on the observed 
properties of stellar populations is that their structural properties 
vary considerably with metallicity and abundance. This has been 
pointed out by |Bovy et aL] ( |2012c| ), who found that when they 
sliced the SEGUE G-dwarf sample into "mono- abundance" bins, 
the structural parameters varied smoothly from short and thick at 
the metal-poor a-rich end, to long and thin at the metal-rich a- 
poor end. We show a similar dissection of our model disc stars in 
Fig. [7] and find a qualitatively similar result, i.e. the scale height 
decreases smoothly from upper-left to lower-right and the thinner 
populations tend to have longer radial scale lengths. 

To obtain the structural parameters, we parametrize the stellar 
distribution as 

p(r,z) oc f(r)f(z), 
where the radial and vertical density functions are 



f(r) exp 




and 

/(,) = sech 2 (-J!-), 

respectively. We fit these functions to the simulated particle density 
distribution for the maximum likelihood values of h r and h z using 
a procedure similar to the one described in Sect. |3.2| We use stars 
within a range of 4 < R < 9 kpc and \z\ < 3 kpc to obtain the fits. 
By requiring at least 100 particles per each plotted cell, we obtain 
lcr uncertainties that are everywhere less than 10%. As in Sect. |3.2| 
we experimented also with an exponential function for the vertical 
distribution, and found that we can obtain smaller errors using the 
sech 2 form. If we fit the vertical distributions of stellar populations 



in our model with a single exponential, we obtain scale heights that 
are somewhat higher than those shown in Fig. [7] similar to what we 
found in Sect.EU 



The locus of "saturated" scale-lengths in Figure [7] below 
[O/Fe] rsj —0.05 are all young stars. Their scale-lengths are long 
because they are forming out of the cold gas which is also very ex- 
tended. A similar saturation is seen in |Bovy et al.| ( |2012c] ). In our 
model, the sub-solar part of this locus originates in the outer disc, 
and |Bovy et al.| { [2012c| ) similarly find that the mean radii for these 
stars are beyond the solar neighborhood. These are young, 0-2 Gyr 
old stars that formed in the outer disc and scattered into the solar 
neighborhood (see also Haywood 2008). They did not migrate via 
the corotation mechanism, because their tangential velocities ex- 
ceed the circular velocity by ~ 20 km/s. Therefore, they heated 
and although one would naively expect their negative AR to lead 
to a thinner scale height, their scale heights are slightly higher than 
the local young population by virtue of this heating. We can identify 
this population in the lower right panel of Fig. [2] as the 2 Gyr-old 
stars that migrate inward by ~ 1 — 2 kpc and find that their scale 
height h z ~ 0.3 kpc. 



An interesting discrepancy with the |Bovy"eT al. (2012c) dis- 
tribution is seen at the extreme metal-rich, a-poor corner ([O/Fe] 
< —0.1 and [Fe/H] > 0.1). We find that those stars have a short 
scale-length much like the older, more a-rich population. In Bovy 
|et al.| < [2012c| ) no such population exists. In our model, this is 
a unique population: comparing their properties with Fig. 14 of 
Loebman et al. ( 2011), we see that these are intermediate- age stars 
(5-6 Gyr old) that come from the same part of the disc as some of 
the oldest stars in the solar neighborhood (Rform ~ 2 kpc). These 
are the prototypical "extreme" migrators that were able to migrate 
far because they never heated considerably. Hence, their velocity 
dispersion is ~ 10 km/s lower than other stars of similar age and 
because of the lower heating they remain thinner than the average 
population of the same age. We speculate that should such stars ex- 
ist in the MW, they may be absent from the SEGUE sample due to 
their thin distribution and low numbers. However, should they be 
observed in future surveys, they would be an indication of signifi- 
cant radial mixing in the MW disc. 
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Figure 7. The distribution of stars in the [O/Fe] vs. [Fe/H] plane, coloured by scale length h r on the left and scale height h z on the right. The logarithmically- 
spaced contours indicate mass density from 10- 10 4 particles per bin. The stars lie between 4 < R [kpc] < 9 and \z\ < 3 [kpc]. 



6 CONCLUSIONS 

We have demonstrated that an outward migrating population of 
stars in a galactic disc always thickens. At the same time, the verti- 
cal velocity dispersion of such a population decreases, but we find 
that despite this decrease the old stellar populations arriving from 
the interior at the solar neighborhood match the kinematics of the 
MW thick disc. Further, we find that radial migration and internal 
heating thicken coeval stellar populations by comparable amounts. 
The thickening due to radial migration alone is well- approximated 
by a simple analytic treatment that assumes a conservation of ver- 
tical action. Importantly, we find that while radial migration does 
cause an increase in scale height with radius, the flaring that results 
is minor. If radial migration can contribute to a thickened compo- 
nent of a galactic disc then it has important consequences in the 
broader context of disc galaxy formation since most other mecha- 
nisms proposed to form a thick disc component involve the cosmo- 
logical environment. 

Our simulation recovers the qualitative structural stellar pop- 
ulation trends observed in the MW SEGUE data. In particular, we 
find a similar dependence of scale height and scale length on oxy- 
gen abundance and metalicity to those found in Bovy et al. ( 2 012c| ). 
However, our model produces a thick disc that is somewhat too thin 
compared to the MW. This is to be expected, since the disc we con- 
sidered in this work is built up entirely from internal mechanisms 
in the absence of a cosmological environment. Additional pertur- 
bations in a more realistic setting would thicken the disc further. 

Finally, we identify an interesting sub-population in the low- 
[O/Fe] high- [Fe/H] corner of the abundance plane (Fig. [7}, which 
appears to be absent from the |Bovy et al.| ( |20i2c] > distribution, per- 
haps due to the lack of low-latitude coverage in the SEGUE data. 
These stars are the prototypical extreme-migrators having migrated 
from ~ 2 kpc in a few Gyr. 
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